library(tidyverse)
library(terra)
library(sf)
library(leaflet)
library(whitebox)
library(scico)Data for testing IHS
Calculating covariates
May as well start with the same set as the source papers - slope, curvatures, TWI ¯\_(ツ)_/¯
Using whitebox here as it’s geomorphometry algorithms apparently handle outliers/noise better in DEMs with projected coordinate systems vs the algorithms implemented in terra::terrain().
elev_dir <- file.path('data_spatial', 'elevation')
if(!dir.exists(file.path('data_spatial', 'terrain_morphometry'))) {
dir.create(file.path('data_spatial', 'terrain_morphometry'))
}
terr_dir <- file.path('data_spatial', 'terrain_morphometry')
wbt_slope(
dem = file.path(elev_dir, 'Lidar_FPS_1m.tif'),
output = file.path(terr_dir, 'Slope_1m.tif'),
units = "degrees",
wd = getwd(),
compress_rasters = TRUE)
wbt_plan_curvature(
dem = file.path(elev_dir, 'Lidar_FPS_1m.tif'),
output = file.path(terr_dir, 'PlanC_1m.tif'),
log = FALSE,
wd = getwd(),
compress_rasters = TRUE)
wbt_profile_curvature(
dem = file.path(elev_dir, 'Lidar_FPS_1m.tif'),
output = file.path(terr_dir, 'ProfC_1m.tif'),
log = FALSE,
wd = getwd(),
compress_rasters = TRUE
)
slope <- rast(file.path(terr_dir, 'Slope_1m.tif'))
planc <- rast(file.path(terr_dir, 'PlanC_1m.tif'))
profc <- rast(file.path(terr_dir, 'ProfC_1m.tif'))TWI is a little more complicated as the DEM needs to be hydrologically corrected first.
# breach depressions (less interferey than filling sinks)
wbt_breach_depressions_least_cost(
dem = file.path(elev_dir, 'Lidar_FPS_1m.tif'),
output = file.path(elev_dir, 'Lidar_FPS_BD_1m.tif'),
dist = 250,
max_cost = NULL,
min_dist = TRUE,
flat_increment = NULL,
fill = TRUE,
wd = getwd(),
compress_rasters = FALSE
)
dem <- rast(file.path(elev_dir, 'Lidar_FPS_1m.tif'))
dem_bd <- rast(file.path(elev_dir, 'Lidar_FPS_BD_1m.tif'))
# check diff
dem_diff <- dem_bd - dem
dem_diff[dem_diff == 0] <- NADifferences after breaching range between -2.21 m and 0.2 m and are largely confined to single-pixel channels in open water, so that’s fine.
# for TWI, slope needs to also come from the breached DEM
# for consistency
wbt_slope(
dem = file.path(elev_dir, 'Lidar_FPS_BD_1m.tif'),
output = file.path(terr_dir, 'Slope_BD_1m.tif'),
units = "degrees",
wd = getwd(),
compress_rasters = TRUE)
wbt_d_inf_flow_accumulation(
input = file.path(elev_dir, 'Lidar_FPS_BD_1m.tif'),
output = file.path(terr_dir, 'Dinf_SCA_1m.tif'),
out_type = "sca",
threshold = NULL,
log = FALSE,
clip = FALSE,
pntr = FALSE,
wd = getwd(),
compress_rasters = TRUE)
wbt_wetness_index(
sca = file.path(terr_dir, 'Dinf_SCA_1m.tif'),
slope = file.path(terr_dir, 'Slope_BD_1m.tif'),
output = file.path(terr_dir, 'TWI_1m.tif'),
wd = getwd(),
compress_rasters = TRUE
)
twi <- rast(file.path(terr_dir, 'TWI_1m.tif'))Not so sure I like this output, but it gives us something to cluster…